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Owing to the recent, rapid development of computer technology, the resolution of atmospheric numerical models has 
increased substantially. With the use of next-generation supercomputers, atmospheric simulations using horizontal grid 
intervals of 0(100) m or less will gain popularity. At such high resolution more of the steep gradients in mountainous 
terrain will be resolved, which may result in large truncation errors in those models using terrain-following coordinates. 
In this study, a new 3D Cartesian coordinate non-hydrostatic atmospheric model is developed. A cut-cell representation 
of topography based on finite-volume discretization is combined with a cell-merging approach, in which small cut-cells 
are merged with neighboring cells either vertically or horizontally. In addition, a block-structured mesh-refinement 
technique is introduced to achieve a variable resolution on the model grid with the finest resolution occurring close to 
the terrain surface. The model successfully reproduces a flow over a 3D bell-shaped hill that shows a good agreement 
with the flow predicted by the linear theory. The ability of the model to simulate flows over steep terrain is demonstrated 
using a hemisphere-shaped hill where the maximum slope angle is resolved at 71°. The advantage of a locally refined 
grid around a 3D hill, with cut-cells at the terrain surface, is also demonstrated using the hemisphere-shaped hill. The 
model reproduces smooth mountain waves propagating over varying grid resolution without introducing large errors 
associated with the change of mesh resolution. At the same time, the model shows a good scalability on a locally 
refined grid with the use of OpenMP. 
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I. INTRODUCTION 

One of the pressing concerns of next-generation high- 
resolution atmospheric modeling is the accurate treatment of 
terrain. The continuous increase in computer power and the 
associated increase in model resolution has resulted in the 
resolution of steeper and more complex features in the ter¬ 
rain. These variations in the bottom surface of the atmosphere 
not only have a significant influence on the local dynamics 
near the surface but can also affect the global circulations 
(McFarlane 1987). Although the resolutions of 1-20 km are 
used in today’s operational models, higher-resolution simu¬ 
lations at horizontal grid intervals of 0(100) m or less will 
gain popularity with the use of next-generation supercomput¬ 
ers (e.g., Miyamoto et al. 2013). Therefore it becomes more 
important for the future high-resolution models to implement 
a robust method of representing topography for steep gradients 
and complex geometries. 

For many years, the common choice for the represen¬ 
tation of topography in atmospheric models has been the 
terrain-following vertical coordinates based either on pres¬ 
sure (e.g., Phillips 1957, Simmons and Burridge 1981), or on 
height (e.g., Gal-Chen and Somerville 1975). In the terrain¬ 
following coordinates, the vertical model levels follow the 
shape of the terrain at the bottom and gradually revert to 
horizontal surfaces with increasing height above the surface. 
The main advantage of this approach is that the imposition 
of the lower boundary condition is straightforward for arbi¬ 
trary topography. In addition, the terrain-following coordi- 
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nates are suitable for coupling with boundary layer param- 
eterizations because a high near-ground resolution is easily 
achieved by increasing the number of model levels near the 
bottom boundary. Although the terrain-following coordinates 
have proven effective for a wide range of applications, large 
truncation errors may arise in computing the horizontal gra¬ 
dients, particularly in the presence of steep terrain (Satomura 
1989, Thompson et al. 1985). It is recognized that the most 
critical error lies in the discretization of the horizontal pres¬ 
sure gradient term, which can induce spurious circulations 
over mountainous topography (Janjic 1989) and even numeri¬ 
cal instability if the mountains are steep enough (Zangl 2012). 

There are many ongoing efforts to alleviate the disadvan¬ 
tages of the terrain-following coordinates. A substantial im¬ 
provement has been made by the specification of the vertical 
coordinate which removes the influence of the topography as 
fast as possible with height (Klemp 2011, Leuenberger et al. 
2010, Sch^etal. 2002). This successfully reduced the er¬ 
rors at upper levels where the coordinates are much smoother 
than classical hybrid coordinates. Another improvement lies 
in a better treatment of the horizontal gradient of pressure 
(Klemp 2011, Zangl 2012) as well as of diffusion (Zangl 
2002). While these studies have improved the accuracy of 
the terrain-following coordinates substantially, a conclusion 
has not been reached on whether the terrain-following mod¬ 
els would be accurate enough for future generation high- 
resolution models. 

This study examines methods based on the use of Cartesian 
coordinates as an alternative means of representing topogra¬ 
phy. In these methods, the terrain is directly incorporated into 
a regular rectangular grid without using a coordinate transfor¬ 
mation. Hereafter in this study, a regular rectangular grid sys- 
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tem is referred to as a ‘Cartesian grid’ and a representation 
method of topography based on a Cartesian grid is referred 
to as a ‘Cartesian-grid method’. Since the model levels are 
kept horizontal throughout the domain. Cartesian-grid meth¬ 
ods resolve the imbalances which occur on a terrain-following 
grid in the discretization of the horizontal gradients. However 
the imposition of the lower boundary condition can be com¬ 
plicated in Cartesian-grid methods because the terrain surface 
does not normally coincide with the grid lines. 

The step-mountain method is a Cartesian-grid method 
with a straightforward imposition of the lower boundary 
condition by approximating the terrain surface as a piece- 
wise constant function along the grid lines (e.g., Bryan 
1969, Mesinger et al. 1988). Because of the lack of ac¬ 
curacy in the stepwise approximation of topography, how¬ 
ever, this method introduces serious errors at step-corners 
into flow patterns (Callus and Klemp 2000) and therefore 
turned out to be ill suited for high-resolution simulations 
over mountains (Schar et al. 2002, Zangl 2003). The 
partial-cell method provides another stepwise approximation 
where the heights of the steps are adjusted to those of 
the topography (e.g., Maier-Reimer and Mikolajewicz 1992, 
Semtner and Mintz 1977). Adcroftetal. (1997) showed that 
the partial-cell method successfully reduced the errors associ¬ 
ated with the stepwise boundary over gently-sloping terrain. 
However the errors can still be large at the step-corners over 
steep terrain where the topographic height varies substantially 
in a horizontal grid length (Yamazaki and Satomura 2010). 

A smoother and more precise representation of the terrain is 
achieved by allowing linear variation of the boundary within 
a cell, resulting in various shapes of cells that are cut by the 
terrain surface. Finite-volume discretization of the governing 
equations assures conservation of model variables on those 
irregularly-shaped cells as well as on regular uncut cells. This 
finite-volume based Cartesian-grid method is refered to as the 
cut-cell (or shaved-cell) method, and it is the approach ex¬ 
plored in this paper. 

The cut-cell method has been most popular in the held 
of computational fluid dynamics for simulating flows with 
complex geometry (e.g., Pemberetal. 1995, Quirk 1994, 
Udaykumar et al. 1996, Ye et al. 1999). After initial imple¬ 
mentation in an ocean model (Adcroft et al. 1997), the cut-cell 
method has been examined for accuracy and robustness in at¬ 
mospheric models over the past dozen years (e.g.. Good et al. 
2013, Jebens et al. 2011, Klein et al. 2009, Locketal. 
2012, Steppeler et al. 2002’ 2006, Walko and Avissar 2008, 
Yamazaki and Satomura 2008’ 2010). Many applications to 
well-known idealized flows in these studies suggest that the 
cut-cell method does not suffer from the problems reported by 
Gallus and Klemp (2000) in the step-mountain method and 
can reproduce smooth mountain waves over the terrain. 

Comparisons of 2D flow results using the cut-cell method 
to results from terrain-following models have been made in 
various studies. For example, Yamazaki and Satomura (2008’ 
2010) demonstrated that the cut-cell method successfully elim¬ 
inated the spurious vertical velocity modes that occurred 
in the vicinity of steep slopes in a terrain-following model. 
Good et al. (2013) compared the errors in the flow aloft, where 
grids are fully rectangular in cut-cell models but are distorted 
in terrain-following models because of the influence of the un¬ 


derlying topography. They showed that the errors associated 
with a terrain-following grid are reduced when the cut-cell 
method is used. They also examined the robustness of the 
cut-cell method for steep gradients and demonstrated that it 
can produce stable results for flows over bell-shaped hills with 
aspect ratios of the height to the half-width up to 10. 

The main problem associated with the use of the cut-cell 
method is the generation of arbitrarily small cut-cells near the 
boundary. Such cells lead to severe stability constraints as a 
result of the Courant-Friedrichs-Lewy (CFL) condition, and 
therefore require very small time steps. Several approaches 
have been introduced to resolve this small-cell problem in at¬ 
mospheric cut-cell models. In a simple approach proposed by 
Steppeler et al. (2002), the computational volumes of cut-cells 
are artificially increased to those of regular uncut cells. This 
method was called ‘thin-walF approximation because, by ex¬ 
panding the volume to the notionally full value but leaving the 
areas untouched, the terrain looks like a collection of infinites¬ 
imally thin-walls (Adcroft 2013). Yamazaki and Satomura 
(2010) use a different approach in which small cut-cells are 
merged with adjacent cells either vertically or horizontally. 
This cell-merging technique makes it possible to extend the 
stability limit maintaining the rigid evaluation of cut-cell vol¬ 
umes and areas, and thus maintaining a sharp representation of 
the terrain surface. In Klein et al. (2009), a dimensional split¬ 
ting technique was used to approximate the fluxes at cut-cell 
interfaces, thereby allowing the use of a full-time step defined 
by the regular grid. The use of implicit schemes, introduced in 
an atmospheric model by Jebens et al. (2011), is another way 
to stabilize small cells. 

Another disadvantage of the cut-cell method compared to 
the terrain-following approach is that to obtain a high vertical 
resolution at boundary layers over a wide range of topographic 
height can be expensive (Walko and Avissar 2008). In terrain¬ 
following models, a high vertical resolution near the terrain 
surface is easily achieved over all topographic heights by in¬ 
creasing the number of model levels near the bottom bound¬ 
ary. However, in a Cartesian-grid based model, a substantial 
number of additional model levels would be required to cover 
all topographic heights. Around steep slopes, horizontal grid 
intervals must be closely spaced as well as vertical grid in¬ 
tervals to achieve high near-ground resolution on a Cartesian 
grid. A block-structured mesh refinement approach proposed 
in Yamazaki and Satomura (2012) is one way to achieve a lo¬ 
cally refined Cartesian grid with high computational efficiency. 
They demonstrated that a flow over a 2D semicircular hill was 
successfully reproduced on a grid locally refined around the 
hill with the use of cut-cells at the boundary. 

Compared to the number of 2D studies of the cut-cell rep¬ 
resentation of topography, applications of cut-cells to 3D at¬ 
mospheric model studies are relatively scarce. Steppeler et al. 
(2006) investigated the impact and potential use of cut-cells in 
a 3D forecasting model. They demonstrated that precipitation 
scores and RMSE of the temperature for 1-day forecasts using 
a cut-cell model were improved compared to the results from 
the terrain-following version of the model. In the extended 
forecasts of 5 days, the improvements became more substan¬ 
tial (Steppeler et al. 201T 2013). 

Lock et al. (2012) focused on the capability of the cut-cell 
method for 3D idealized flows with a very steep slope. They 
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showed that a potential flow over a 3D bell-shaped hill was 
successfully reproduced using cut-cells where the maximum 
slope angle of the hill was as steep as 74° from the horizon¬ 
tal. However no method to address the small cell problem was 
adopted in their cut-cell model, leaving stability and efficiency 
issues for practical applications. 

In this study, we propose a new 3D atmospheric model us¬ 
ing the cut-cell method for high-resolution simulations over 
steep topography. To avoid the severe stability constraints 
from small cells, the cell-merging technique developed by 
Yamazaki and Satomura (2010), hereafter YSIO, is extended 
to 3D and implemented in the model. In addition, a 3D version 
of the block-structured mesh refinement approach proposed in 
Yamazaki and Satomura (2012), hereafter YS12, is introduced 
to achieve high resolution near the terrain surface and also easy 
parallelization. Section 2 describes how the 2D methods of 
YSIO and YS12 are converted into 3D. Then we demonstrate 
the performance of the model in section 3 through numerical 
experiments on mountain waves. To investigate the capability 
of the proposed 3D cut-cell method for a wide range of slope 
angles, the results of flow over a 3D bell-shaped hill and a 
hemisphere-shaped hill are presented. In addition, the perfor¬ 
mance of the model on a locally refined grid in combination 
with the use of cut-cells near the boundary is examined using 
the hemisphere-shaped hill. 


gas constant, and po is a reference pressure of 10^ Pa. The 
terms Du, Dy, Dyj and Dg in Eqs (l)-(4) represent source 
terms due to mixing and diffusion. In this study, turbu¬ 
lent parameterization terms based on turbulent kinetic energy 
(Klemp and Wilhelmson 1978) are used. In addition, a fourth- 
order artificial diffusion term is introduced in the horizontal 
and vertical directions to suppress the numerical noise. Finally 
the system of equations is closed by the following equation of 
state for an ideal gas: 


6 = 


P 

pR 



(9) 


In the view of conservation characteristics, flux-form equa¬ 
tions are well suited to the finite-volume discretization that is 
an essential part of the cut-cell method. Satomura and Akiba 
(2003) demonstrated that the equations achieve mass conser¬ 
vation by simulating heat-island circulation. In addition, they 
designed the equations to avoid cancellation errors stemming 
from subtracting the hydrostatic variable (p or p) from the 
nearly hydrostatic total variable (p or p). This may occur in 
other flux-form equations (e.g., Klemp et al. 2007, Saito et al. 
2001). Specifically, Satomura and Akiba (2003) accomplish it 
by directly predicting the perturbations of the variables (p' or 
P')- 


II. MODEL DESCRIPTION 


B. Block-structured grid 


A. Governing equations 


The model solves fully compressible quasi-flux-form equa¬ 
tions developed by Akiba (2002) and Satomura and Akiba 
(2003), given by the following conservation equations for mo¬ 
mentum, potential temperature and mass based on the Carte¬ 
sian coordinates: 

Qn' 

= -V • (puu) “ > (1) 

dp' 

= -V • (pvu) - ^ + Dy , (2) 

oy 

9 75 ^ 

= -V • (pwu) - — - p'g + Du, , (3) 

CvPO \PoJ 

= -V • (pu) , (5) 

where u = {u, v, w) is the velocity vector, 6 is the potential 
temperature, g is the acceleration due to gravity, p and p are the 
total pressure and total density, respectively, and the prime in¬ 
dicates perturbations from the hydrostatically balanced state: 


dpu 

dt 

dpv 

dt 

dpw 

dt 

dp' 

dt 

V 

dt 


P — P{z) + P[x,y,z,t)-< 
P = P{z) + P'(x,y,z,tp 


dp 

Vz 


= -P9- 


( 6 ) 

(7) 

( 8 ) 


In Eq. (4), Cp and Cy denote the specific heats at con¬ 
stant pressure and constant volume, respectively, R is the 


To achieve variable resolution in the Cartesian-coordinate 
system, the generation of the model grid begins with the 
creation of a block-structured grid. Here we use the Con¬ 
served Building-Cube Method (CBCM) that is originally pro¬ 
posed in 2D by YS12. The method is based on the Bulding- 
Cube Method (BCM) developed by Nakahashi (2003) that has 
been used for the problems of computational fluid dynamics 
(e.g., Kim et al. 2007, Nakahashi et al. 2006). CBCM and 
BCM share the two-tiered data structure of a generated block- 
structured grid: tree-structured sub-domains called ‘cubes’ 
and array-structured uniform mesh called ‘cells’ in each cube. 

Figure 1 illustrates a 2D example of the grid generation pro¬ 
cess of CBCM. First the model domain is divided into coarse 
equally-spaced cubes (Figure 1(a)). Next, cubes that are in the 
vicinity of the boundary are divided into four cubes (or eight 
in 3D). By repeating this refinement process, the domain is 
divided into a number of cubes, where the size of a cube be¬ 
comes small closer to the terrain surface (Figure 1(b)). Here 
the size differences between cubes are adjusted to guarantee 
a uniform 2 : 1 mesh resolution at fine-coarse cube bound¬ 
aries in horizontal, vertical and diagonal directions. Note that 
the cubes that are located completely inside of the topography 
are removed and are not used for the computation. Finally, a 
Cartesian grid of equal spacing and equal number of cells in 
each cube is generated (Figure 1(c)). For example, each cube 
in Figure 1(c) has 4^ cells (or 4^ in 3D) regardless of the size 
of the cube. Note that the generation of cut-cells in the finest 
cubes is described in the following subsection IIC. 

The local grid interval is determined as 

H{1) = 2-^H, 


( 10 ) 
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(a) (b) (c) 

FIG. 1. Schematics of the generation of a block-structured Cartesian mesh: (a) and (b) show the cubes before and after the refinement process, 
respectively, and (c) shows the cells as well as the cube boundaries. Thick and thin lines represent cube boundaries and cell boundaries, 
respectively. Shaded regions describe the topography. 




where the refinement level I ranges from 0 (coarsest) to Imax 
(finest) and the grid interval at the unrefined cube H = H{0). 
In this study, the same grid interval H is used in x, y and z 
directions for simplicity, though it can easily be extended to 
different spacing in each direction. When setting up a block- 
structured mesh, we choose the value of H so that the field far 
away from the topography is resolved at this resolution. Then 
we keep refining the cubes near the terrain surface until we 
obtain as high a resolution as we wish to impose close to the 
surface. For example, when we use H = I km, we can get a 
local grid interval of 15.625 m at the refinement level 6. 

A block-structured grid constructed from cubes and cells 
provides several attractive features for Cartesian-grid models. 
First, the tree-based data structure of cubes provides a locally 
refined grid around arbitrary topography by refining the size of 
the cubes near the terrain surface. The use of a uniform Carte¬ 
sian mesh in each cube, at the same time, allows the direct 
use of any existing code based on a Cartesian grid by treat¬ 
ing each cube as an independent computational domain. Fur¬ 
thermore, the same number of cells among all cubes makes 
the method suitable for parallel computing. Because the load 
balance of each cube is equivalent, high parallel efficiency 
can be achieved by simply distributing an equal number of 
cubes for each processor. Both CBCM and BCM showed good 
speedups with straightforward OpenMP-based parallelization 
(Kim et al. 2007, Yamazaki and Satomura 2012). 

CBCM employs a subcycling time integration that allows 
the use of a larger time step at coarse cubes than that used for 
smaller cubes. For cubes at the refinement level I, the time 
step is chosen as 

At{l) = 2-^At, (11) 

where the global time step At = At{0). The equations are 
integrated from the cubes at the finest level. After the cells in 
the cubes at level I are advanced in two time steps of At{l), 
those at level Z — 1 are advanced in one time step of At {I — 
1), followed by the information exchange between the levels. 
In our model, this subcycling integration is incorporated into 
the leap-frog time-stepping scheme and used along with the 
Robert-Asselin filter (Asselin 1972, Robert 1966). 


How information is transferred between adjacent cubes 
through ghost cells that are added beyond the boundary of 
each cube. In this study, four ghost cells are added as shown 
in Figure 2. Note that all ghost regions are at the same reso¬ 
lution as the inner domain of the cube. The information be¬ 
tween the same-size cubes can be exchanged in a straightfor¬ 
ward way because of the exact overlapping of the cells (Figure 
2(a)). Some interpolation methods are required at fine-coarse 
cube boundaries to exchange the boundary values at different 
resolutions (Figure 2(b)). In this study, the values of ghost 
cells of a coarse cube are assigned by using a simple four- or 
eight-point average in 2D or 3D, respectively, of fine-cell val¬ 
ues at the corresponding location in the fine cube. The ghost¬ 
cell values of a fine cube, on the other hand, are interpolated 
by assigning the same coarse-cell value to the corresponding 
four or eight ghost cells of the fine cube in 2D or 3D, respec¬ 
tively. Though this interpolation procedure has been demon¬ 
strated in YS12 to conserve the global second-order accuracy 
of the CBCM, a higher-order interpolation method can be used 
to improve the accuracy at fine-coarse cube boundaries (e.g., 
Jablonowski et al. 2006). 

Without time interpolation, some time inconsistency may 
occur during the subcycling integration. For example, in a 2D 
case of Figure 3, the coarse-cell value at (c) is assigned to the 
ghost cells of the fine cube (d) to (g) and can be used twice 
in a row without updating. This causes computational modes 
to the fine-coarse cube boundary due to the time suspension 
at the ghost region. CBCM avoids this problem by integrating 
equations on some of the ghost cells of the fine cube and us¬ 
ing the updated values for the subcycling steps. The number 
of ghost cells on which the equations are integrated is chosen 
to satisfy the numerical stencils at the boundary cells inside 
the fine cube, such as the cells (a) and (b) in Figure 3. In our 
model, we use a 5-point stencil in the horizontal and vertical 
directions to calculate fourth-order artificial diffusion terms. 
Therefore we integrate equations on two more cells beyond 
the boundary to satisfy the numerical stencils of the boundary 
cells with updated values. In case of Figure 3, the ghost-cell 
values at (d) to (g) are updated as well as the values at (a) and 




























































































































































5 


(a) 



FIG. 2. Information transfer between the adjacent cubes of (a) the 
same size and (b) different sizes. Thick lines represent cube bound¬ 
aries, and shaded regions describe the overlapped regions of the ghost 
cells. Solid arrows in (a) describe information transfer between the 
same-size cubes. Thick solid and dotted arrows in (b) describe infor¬ 
mation transfer from a coarse cube to a fine cube, and transfer from a 
fine cube to a coarse cube, respectively. In this figure a cube has 8 x 
8 cells with 4 ghost cells beyond the boundary. 


(b) during the subcycling integration. To be able to integrate 
equations on two more cells beyond the boundary, four ghost 
cells in total are used in this study beyond each side of the cube 
boundary. This approach prevents the computational modes 
from contaminating the results inside the cubes by simply in¬ 
tegrating equations on some extra cells beyond the boundary 
of each cube. On the other hand, it demands higher computa¬ 
tional costs compared to a time interpolation scheme due to a 
relatively large number of ghost cells. 

Another characteristic of CBCM is that it ensures global 
mass-conservation on a locally rehned mesh. To achieve con¬ 
servation with a subcycling time-stepping scheme, we must en¬ 
sure that, at each hne-coarse grid interface, the numerical flux 
on the coarse grid equals the flux on the fine grid accumulated 
during the subcycling steps. CBCM facilitates this process by 
introducing the cube-boundary flux at fine-coarse cube bound¬ 
aries. For example, in the case of Figure 3, the cube-boundary 
flux Fab is defined at the same location as the coarse-cell flux 
fc- During the subcycling integration of the fine cube, the fine- 
cell fluxes fa and fb are accumulated and stored in Fab- Then, 
in the integration of the coarse cube. Fab overrides the coarse 
cell flux fc and is used to update the coarse-cell value at (c). 
In the 3D method, four fine-cell fluxes are accumulated and 
stored in each corresponding cube-boundary flux. This flux¬ 
matching algorithm ensures mass-conservation on the condi¬ 
tion that the density is defined at the cell-centers, as demon¬ 
strated in YS12. 


fine-coarse cube boundary 




FIG. 3. Computational cells and fluxes at a fine-coarse cube bound¬ 
ary. Thick and thin lines represent cube and cell boundaries, respec¬ 
tively. Dotted lines represent the boundaries of ghost-cells of the fine 
cube. Filled and open circles describe computational nodes of the 
cells inside the cube boundary and those of the ghost-cells, respec¬ 
tively. Thin and thick arrows indicate the fluxes at the cell and cube 
boundaries, respectively. 


C. Cut-cell configuration 


Following the generation of a block-structured grid in the 
previous subsection, cut-cells are generated near the terrain 
surface. Because the grid is refined around the topography, 
the procedure of cut-cell generation is only required in the 
finest cubes. Following Steppeler et al. (2006) and Lock et al. 
( 2012 ), the topographic boundary is represented by piecewise 
bilinear surfaces that are continuous at the boundaries of grid 
columns. First the terrain heights are specified at the four cor¬ 
ners of the grid columns in each cube, then a unique surface 
for each grid column is defined by using a bilinear function of 
height with respect to the horizontal position {x, y) as 

h(x, y) = mix + m 2 Xy + m^y + c, ( 12 ) 

where mi, m 2 , m 3 and c are constants. The function gives a 
linear spline at any vertical cross-sections in the x or y direc¬ 
tion, as shown in Figure 4. By treating the corner O at (z —1/2, 
j — 1 / 2 ) as the origin, the four heights at the corners (z ± 1 / 2 , 
j ± 1 / 2 ) determine a bilinear surface at the grid column (z, j) 
through the coefficients 


mi,, = 


m 2 ,, = 


"i3,, = 


~ 2-^2 2'J 2 

Ax 

~ ^i-y+h + ^i+y+h ~ ^^+y-h 


Ax Ay 


* 2 -'^ 2 2 ^ 2 


Ay 


Cij — 


(13) 

L 

1 

(14) 

(15) 

(16) 


where Ax and Ay indicate the grid intervals in the x and y 
directions, respectively, which are equal to H{lmax) in this 
study. 

Based on the bilinear representation of the terrain surface, 
the volumes and areas of the cut-cells are computed. To en¬ 
able easy computation of the 3D cut-cell parameters, we use 
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FIG. 4. Bilinear representation of the terrain surface at the grid col¬ 
umn {i, j). The corner O acts as the origin of the bilinear function at 
the column. Thick lines indicate the boundaries at the vertical cross 
sections through the center of the grid column, which determine the 
gradients of the surface. 


dh/dy 



□ 


\~D 


vertical merging 


horizontal merging 
(±x direction) 

horizontal merging 
(±y direction) 


FIG. 5. The direction of cell-merging depending on the horizontal 
gradients of the bilinear surface. The shaded region (a) describes the 
range of the gradients where a small cell is merged with each upper 
cell. Hatched regions describe the range of the gradients where small 
cells are merged with each adjacent cell in the (b) +x and (c) —x 
directions, respectively, and dotted regions describe the same in the 
(d) +y and (e) —y directions, respectively. 


the approach of Lock et al. (2012), which divides a cell into 
narrow rectangular columns and approximates the cut-cell vol¬ 
ume by the sum of those volumes. The areas of the cut-cells 
are also computed based on the same approximation. Detailed 
explanation of this computation is found in Lock et al. (2012); 
see their Appendix B. 

Cell merging is an important feature in the model described 
in this study. In this method, cut-cells whose center is under¬ 
ground or whose volume is smaller than half the volume of a 
regular cell are merged with an adjacent cell either vertically 
or horizontally. The direction of cell-merging is determined by 
the horizontal gradients of bilinear surfaces that are evaluated 
at the center of each grid column. In the case of Figure 4, the 
gradients of the surface in the x and y directions are evaluated 
at the grid column (i, j) as 





9a; y 

). - 

2Ax 


(17) 

9/i^ 



92/y 

) - 

to 

00 


respectively. When both \{'dh/dx)ij\ and \{dh/dy)ij\ are less 
than or equal to 1, small cut-cells at the column (i, j) are 
merged vertically with each upper cell. Otherwise, they are 
merged with an adjacent cell in one of four horizontal direc¬ 
tions, -\-x, —X, +y or —y, determined by the values of the 
gradients (Figure 5). Note that we switch the direction of cell¬ 
merging between the vertical direction and the horizontal di¬ 
rections at the gradients of ±1. Following the result of YSIO 
that the vertical and horizontal merging of cells gave consis¬ 
tent results of flow over a pyramidal mountain, we use the 
vertical merging at those gradients at our own choice. Simi¬ 
larly, when \{dh/dx)ij \ and \{dh/dy)ij \ are larger than 1 and 
\{dh/dx)ij\ = \{dh/dy)ij\, we use the horizontal merging 


in ±x direction rather than that in ±?/ direction. This cell¬ 
merging procedure can be described in a form of the pseu¬ 
docode shown in Figure 6. 

After the cell-merging procedure, the computational vol¬ 
umes of the cells become larger than half of a regular cell, 
therefore allowing the use of up to half of the full-time step 
defined by the regular cell. In this study we assume that the 
topography is well resolved on the grid, and are not concerned 
with the topography such as an extremely steep v-shaped val¬ 
ley, where a small cut-cell may not have a large adjacent cell 
to merge with. 

Finally, the model variables are arranged on the cells. Fol¬ 
lowing the 2D method of YSIO, a semi-staggered arrange¬ 
ment of variables is used in this study: scalar variables (p', 
p' and 0) are arranged on the cell centers, and all the veloc¬ 
ity components {u, v and w) are co-located and arranged on 
the corners of the cells. For an uncut grid cell centered on 
(i,j, k), velocity components are arranged on the eight corners 
at (i± 1/2, j ± 1/2, fe± 1/2) as shown in Figure 7(a). Here the 
square and circles represent the location of the scalar point and 
velocity points on the cell, respectively. For cut-cells, velocity 
components are arranged on the corners of the cells above the 
surface and also on the topographic boundary, as shown for a 
case with and without cell-merging in Figures 7(b) and 7(c), 
respectively. Variables on horizontally merged cells are also 
arranged in the same way. Note that no variables are arranged 
underground. This unique arrangement of variables enables 
a direct evaluation of the boundary velocity at the terrain sur¬ 
face, as in the 2D method of YSIO, thereby simplifying the 
computation of the velocity near the boundary. The details are 
discussed in the following subsection as well as the computa¬ 
tion of the scalar variables on cut-cells. 


D. Spatial discretization 

To solve flows through the irregularly shaped cut-cells, Eqs 
(l)-(5) are discretized in space using a finite-volume approach 
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if {the cell (i, j, k) needs to be merged } then 

if{\{dh/dx)ij \ ^ 1 and \{dh/dy)ij \ ^ l}then 
merge with the cell (i, j,k + 1 ) 

elseif 

if{{dh/dx)ij > 1 and \{dh/dx)ij\ ^ \{dh/dy)ij\} then 
merge with the cell (i — l,j, k) 

elseif {{dh/dx)ij < —1 and \{dh/dx)ij \ ^ \{dh/dy)ij\} then 
merge with the cell {i + l,j, k) 
elseif {{dh/dy)ij > 1 and \{dh/dx)ij \ < \{dh/dy)i^jW then 
merge with the cell {i,j — 1 , k) 

elseif {(dh/dy)ij < —1 and \{dh/dx)ij \ < \{flh/dy)i^j\} then 
merge with the cell (i,j + 1 , k) 

endif 

endif 

endif 


FIG. 6. Pseudocode of the cell-merging algorithm. 




FIG. 7. Variable arrangement on (a) a regular cell, (b) a non-merged cut cell and (c) a vertically merged cut cell. Thin solid and dashed lines 
describe the grid lines, and thick solid and thick-dashed lines describe the boundaries of the cells. Squares and circles represent scalar points 
and velocity points, respectively. Shaded region represents the topographic surface. 


based on the 2D method of YSIO. The method invokes 
Gauss’s divergence theorem, which states that the volume inte¬ 
gral of the vector divergence over a control volume V enclosed 
by the surface S is transformed to a surface integral as 


tegrals over all the external sides of the volume. Introducing 
T/j as the area mean of the component of F normal to the 
side J of the control volume I, the surface integral in Eq. (20) 
becomes 


V-FdV= (bF-ndS, 


(19) 


where F is a flux vector and n is a unit vector pointing along 
the outward normal of the surface S. Assuming that V ■ F is 
constant over a control volume, the vector divergence through 
a discrete control volume I is expressed as 

{V ■F)i = y^ j^F-ndS, (20) 

Si 


where Vj and Sj is the volume and surface of the control vol¬ 
ume I. Then the surface integral is decomposed into the in- 


F-ndS = 


Si 


J ' 


Si,. 


F-ndS = > Fi,jSij. (21) 


J 


where Si,j is the area of the side J. 

Starting with the simplest case, consider a regular uncut cell 
shown in Figure 7(a). Application of Eqs (20) and (21) on the 
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cell (i,j, k) gives 


follows: 


' ^)ijk — \^^x{^xSx')ijk “t“ ^y{^ySy)ijk ^ziFzSz'jijk} 

Vijk 


f ^P'u\ 

1 r 

\ 9^ Ji'fk' 

Vi'j'k' ^ 


'i'k' ^ 


V (^“i+ 1 

Viik V ‘+5 


—zx — y—y ,~xy — z—z 'j 

+ dy{p U V Sy)i'fk' -\-Oz{p U W Szji'fk'l 


^ijk 

+ Py. 








ij — ^ij — ^k 


F, ,-F, 

^3k^\^■k ^3k•\^■k ijk—^ ^3k~ 


Ax 


-I- D* 

' ^ui' j'k' ’ 


(30) 



/ dpv\ 

1 / 

) , ( 22 ) 

V 9^ Ji'fk' 

Vi'j'k' ^ 


-yz-x-x 

i^p V U Oxji'j'k' 


—zx — y—y ,~xy — z—z "I 

+ 5y[p V V Sy)ifj'k' ^0z{p V W Szji'fk'j 


where 


SxWijk = V'i+ijfc - i’i-^jk : (23) 


and we define i5y('!/')ijfc and Sz{'4’)ijk in a similar manner. Here 
Fx, Fy and F^ are the x, y and z components of the flux vector, 
respectively. On a regular cell, the cell volume Vijk and the 
surface areas normal to the x, y and z directions, Sx, Sy and 
Sz, respectively, are computed as 


Vijk 

= AxAyAz, 

(24) 

^±^3k 

= 

(25) 

^ij ± ^ fc 

= AzAx, 

(26) 


= AxAy, 

(27) 


where Az indicates the grid interval in the z direction. Sup¬ 
posing that F is an advective flux F = (pu, where p is the 
scalar quantity, Eq. (22) becomes 


1 C ^^x^^y z ^^y^^zx 

' F)ijk = — s ^x{,4* tr Sx)ijk “f ^y{4^ Sy^ijk 
Viik t 


ijk 


(28) 


where 


^Ijk = (V'i-ijfc + V'i+ijfc) / 2 , (29) 


and we define and pj^ji. in a similar manner. When the 
volume and areas of a regular cell (24)-(27) are assigned, Eq. 
(28) reduces to a centered finite-difference expression for V • 

F. 


Syipn^fk' , 

' ^vi' j'k' ’ 


Ay 

1 


(31) 


{6x{py^w^u^Sxhfk> 


/ dpw\ 

+ Sy{r^wyvySy)z^,,k'+Sz{r^w^w^Sz)i'fk'} 

- 2"''''' - Ti^^'k'g + , (32) 

^ _c^ fPMlV^ {Sx(^^7^Sxhk 
9f / iik \ po J Vijk ^ 


+ 6y{p0^ Sy)ijk + SzipO^'w'^^ Sz)ijk 


- D 


i 


dt 


9 ijk ^ ’ 

1 

Viik 


(33) 


ijk ''ijk 

+ Szirw^yszhk}, 


{SxirFy^Sx)ijk + 6y{pyV^^Sy)zjk 


(34) 


where Vi'j'k' is the volume of a velocity cell centered on 
= (i — 1 / 2 , J — ll2,k — 1 / 2 ), which is equal to 
Vijk when the velocity cell is uncut. The terms 79*, 79*, 79^, 
and Dg are discretized forms of the diffusion terms. In this 
study the second-order central difference scheme is used for 
pressure gradient terms and diffusion terms. For cells with¬ 
out necessary neighboring fluid points for the calculation of 
a fourth-order artificial diffusion term, a second-order term is 
introduced instead, where boundary conditions are used at the 
terrain surface. 

We will now describe how the finite-volume approach is 
implemented in a situation where some of the cells are cut by 
the topography. As a semi-staggered arrangement of the scalar 
and velocity variables is used, the scalar and velocity cells are 
at different locations and are treated separately. First, consider 
the discretization of conservation equations on the scalar cell 
centered at point Pq (Figure 8 ). The cell is vertically merged, 
and enclosed by eight velocity points at Uq to U 7 . Though we 
focus here on the case of a vertically merged cell, the spatial 
discretization on horizontally merged cells is handled in the 
same way. 

As shown in the case of an uncut cell, finite-volume dis¬ 
cretization requires the estimation of surface integrals over 
each face of the cell. For example, the surface integral of the 
flux over the right face of the cell Pq is evaluated as 


Using the notations introduced above, the spatially dis¬ 
cretized forms of Eqs (l)-(5) for an uncut cell are given as 


^Uoi23 


F-ndS = Fx,Sx 


(35) 
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FIG. 8. Flux calculation on a scalar cell cut by the topography. Ar¬ 
rows indicate fluxes through the faces normal to the x direction of the 
scalar cell centered at Pq. Hatched regions indicate the cell faces nor¬ 
mal to the X direction. Thin solid and dashed lines describe the grid 
lines, and thick solid and thick-dashed lines describe the boundaries 
of the cells. Squares and circles represent scalar points and velocity 
points, respectively. Shaded region represents the topographic sur¬ 
face. 


where U 0123 denotes the face enclosed by the points Uq, Ui, 
U 2 and U 3 . In the right-hand side, Sx^ represents the area 
of the face U 0123 , and Fx^ is the area mean of the flux over 
the face (Figure 8 ). In the case of the left face of the cell Pq, 
U 4567 , the flux through the face is composed of two fluxes: the 
flux through the boundary with the cell P 2 and that through the 
boundary with the cell P 3 . Thus the surface integral of the flux 
over the face U 4567 can be decomposed as 

I F 11(15= I F-11(15+ I F-11(15 (36) 

-^1^4567 -7174589 -fUeTSg 

= Fx,5x, + Fx,5x„ (37) 

where 5x2 are the areas of the boundary faces U 4589 

and U 6789 , respectively, and Fx 2 and Fx^ are the area means 
of flux over each of those faces, respectively. 

The evaluation of 77 ^^ ^ Fx 2 and Fx^ requires the area mean 
of the normal velocity and the advected scalar quantity over 
each face. In this study, the normal velocity of these fluxes is 
estimated by a linear interpolation among the velocity values 
at the corners of each face. For example, the normal velocity 
of Fxi is computed as 

Ux^ = (Wue + '“ui + U^2 + ^U 3 )/4- (38) 

The normal velocity of Fx 2 and Fx^ on the merged face U 4567 
are given by 

UX2 = + Mua + 77ug )/4, (39) 

UX3 = + Uu, + Uu, + Uug )/4, (40) 


respectively. Note that the calculation of the velocity values in 
the right-hand sides of the Eqs (38)-(40) is described later in 
this subsection. 

The calculation of the area mean of the scalar quantity over 
the faces can be more complicated. In the case of Figure 8 , the 
scalar quantity in Fx 2 is obtained by a simple linear interpola¬ 
tion between neighboring cell centers as 

0.. = (•/'Pg + 0pJ/2. (41) 

The evaluation of the scalar quantity in Fx^, on the other hand, 
is not straightforward because one of the neighboring cell cen¬ 
ters is underground and absorbed to the cell Pq in the cell¬ 
merging procedure. One way to evaluate it with a high-order 
accuracy is to use a multi-dimensional interpolation or extrap¬ 
olation scheme around the boundary. A disadvantage of the 
scheme is that it generally entails a considerable complexity, 
especially in 3D, to locate the appropriate grid points to form 
a multi-dimensional function for various shapes of cut-cells. 
For computational simplicity, here we use a first-order calcu¬ 
lation as used in the 2D study of YS12. In this method, the 
area mean of the scalar quantity over a face is approximated 
by the simple average of cell center values of the cells exchang¬ 
ing fluxes through the face. In the case of Figure 8 , the scalar 
quantity in Fx^ and Fx^ is therefore evaluated as 

=(</'pg +</>P,)/2, (42) 

0X3 = (</'p„ +</>P3)/2, (43) 

respectively. Note that this approximation does not violate 
mass-conservation. As a result, the surface integrals of the 
flux over the left and right faces of the cell Pq are computed 
as 


f F-n (15 = 

4^X2^X‘2^X2 ^X^^X^3x^-j 

(44) 

U 4567 



f F-11(13 = 

4 ^X 1 '^X \ 3xx J 

(45) 

^Uoi23 




respectively. The surface integrals of the flux over the other 
faces of the cell, U 1357 , U 0246 and U 0145 , are also computed 
in the same way, with zero normal flux assumed at the topo¬ 
graphic boundary. Boundary conditions that involve nonzero 
normal fluxes, for example heat conduction across the sur¬ 
face, can be written as effective volume-mean source terms 
(Adcroft et al. 1997). 

Next we consider the discretization of momentum equations 
on the velocity cells. A key problem here is to accurately eval¬ 
uate the pressure gradient at each velocity cell. Because some 
pressure points are underground in a cut-cell model, a veloc¬ 
ity cell cut by the topography may not have all the necessary 
pressure points for the calculation of the pressure gradient. To 
allow the same pressure gradient calculation on cut-cells as 
is used for a regular cell, Walko and Avissar (2008) predicted 
approximate solutions of underground pressure by assuming 
that the cell volumes and cell areas of cut-cells are uniformly 
occupied by small solids. The assumption makes the shape 
of topography indistinct and thus could affect the flow dynam¬ 
ics near the topography. Ye et al. (1999), on the other hand, 
proposed an approach to express the pressure held near the 
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surface in terms of a polynomial interpolating function, and 
evaluate the gradients on cut-cells based on the interpolating 
function. They used pressure values at available neighboring 
cell-centers and also at the surface to construct an interpolat¬ 
ing function for each cut-cell. This method allows the system¬ 
atic evaluation of the pressure gradient at cut-cells of various 
shapes. However the calculation of surface pressure causes dif¬ 
ficulties in atmospheric simulations because the zero-gradient 
boundary condition on the pressure field does not hold true for 
density-stratified flow. In addition, since geophysical flow is 
driven by slight perturbations of the pressure from the hydro¬ 
statically balanced state, extrapolation of the pressure pertur¬ 
bation to the surface could also affect the local flow dynamics 
as well as use of the underground pressure. 

The commonly-used staggered arrangement of variables in 
atmospheric models introduces another complexity into the 
discretization of momentum equations. On a non-staggered 
grid, not only are the velocity components and scalar vari¬ 
ables co-located, but the position and geometry of the asso¬ 
ciated cells are also identical. With a staggered grid, the scalar 
and velocity cells are at different locations and will generally 
have a different shape when they are cut by the topography. 
In general, a cut cell method for a staggered grid must deal 
with this extra complexity (Kirkpatrick et al. 2003). With a 
cell-merging approach, in particular, another formulation of 
the discretized equations for merged velocity cells is normally 
required. 

The current method provides a way around these complica¬ 
tions by using the non-conventional semi-staggered variable 
arrangement described in section IIC. For example, consider 
the velocity points illustrated in Figure 9, which are arranged 
on the same scalar cells shown in Figure 8 . First, we solve the 
momentum equations on the velocity cells which retain the 
regular rectangular shape. Supposing that the pressure points 
neighboring Pq, Pi, P 2 and P 3 are all available, the velocity 
cells with centers indicated by open circles have all the eight 
pressure values available at the corners of the cells. Therefore 
the velocity values on those cells are obtained using Eqs (30)- 
(32). Once the velocity at the open-circled points is predicted, 
the remaining velocity points are either on the topographic sur¬ 
face or on the merged faces of the scalar cells, indicated by 
filled circles and crossed circles, respectively (Figure 9). The 
velocity on these points are evaluated diagnostically. 

The velocity on the topographic surface is calculated by ap¬ 
plying a velocity boundary condition on the surface. When 
the non-slip boundary condition is imposed, all the boundary 
velocities on the surface are set to zero. With the free-slip 
boundary condition, the boundary velocity is calculated so that 
the component of the velocity that is tangential to the surface 
is preserved near the boundary. Given that the unit normal 
direction to the surface at a boundary point Ub is , the 
boundary velocity is therefore calculated as 

u,, = u,, — (u,, • n,, )n,, , (46) 

where indicates a mean velocity near the boundary point 
Ub that is computed from the predicted velocities on the regu¬ 
lar velocity cells. 

For example, consider the computation of the velocity at 
the boundary point U 7 using the free-slip condition. Suppos¬ 
ing that the scalar cell Pq is at the grid column (i, j), the unit 



FIG. 9. Velocity calculation on cut-cells. Squares represent scalar 
points. Open circles represent velocity points on which model solu¬ 
tions are predicted. Filled circles and crossed circles represent diag¬ 
nosed velocity points on the terrain surface and on the merged faces 
of the scalar cells, respectively. The dashed-dotted line and the arrow 
describe the normal line and the unit normal direction to the surface at 
the boundary point U7, respectively. The hatched region indicates the 
nearest plane for the boundary point U 7 defined by four open-circled 
points of the scalar cell P 2 . 


normal direction to the surface at the boundary point U 7 , 



= (na;,ny,n^), 

(47) 

computed 

as. 




(48) 


J{dh/dx)^ 1 ^+{dh/dy)f ^ ^ + l‘ 

^ ^ 2 -> 2 2-J 2 


-idh/dy)^y i 

(49) 

liy 

J{dh/dx)^ 1 ^+{dh/dy)f ^ ^ + l‘ 

V ^ 2 J 2 * 2-^2 


1 

(50) 


/(0/i/9a;)2 1 ,+{dh/dy)l , ^ + l‘ 

V 2 J 2 2 J 2 


where the gradients of the surface at the boundary point U 7 
are calculated at their horizontal grid positions (f — 1 , j — 1 ). 
Following Eqs (17) and (18), the gradients are obtained as 


/ dh\ _ (hij-i + hij) — {hi-ij-i -f hi-ij) 

(51) 

/ dh\ _ + hij) — (hi-ij-i -f hij-i) 

V9yA-i7-i~ 2Z\J/ 

(52) 

Next, to estimate a mean velocity near the boundary point, 
we define the nearest plane to the point by four open-circled 
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points within a cell. The mean velocity near the point is then 
evaluated at the intersection of the nearest plane and the nor¬ 
mal line to the surface. When the normal line at the boundary 
point U 7 is as described in Figure 9, the velocities at the four 
points U 4 , U 5 , Uii and U 12 are used to evaluate the mean 
velocity at the intersection point N as 


= al3 + a(l - /3)u„^ + (1 - a)/3 
+ (l-a)(l-/3K^^, (53) 

where a and jS are the linear interpolation factors dehned as 


a = 


Ax 


/3 = 


I/n 2/ui 

Ay 


(54) 


As a result, the velocity at the boundary point U 7 is computed 
as 


Uu7 = Un - (U^ • . (55) 

The velocities at the other boundary points, such as at U2, U3 
and Ue, are also calculated in the same way. 

Finally, the velocities at crossed-circled points on the 
merged faces are calculated assuming a linear distribution of 
the velocity over each merged face. For example, the velocity 
at Ug and Ug on the merged face U4567 of the scalar cell Pg 
is calculated by a linear interpolation between the velocities at 
the corners of the merged face as 

Uug = + -7), (56) 

Uug = Uu^C + Uu7(l - 0, (57) 

respectively, where the interpolation factors 7 and C, are de¬ 
hned as 


7 = ^53^ 

^U4 

All the velocity values are calculated without estimating the 
surface or the underground pressure with the following three- 
step calculation of the velocity: 1 ) solve the momentum equa¬ 
tions on regular velocity cells, 2 ) calculate the velocity on the 
terrain surface using a boundary condition, and 3) calculate 
the velocity on merged faces by an interpolation of the veloc¬ 
ity obtained in 1) and 2). In addition, there is no need to merge 
velocity cells, and thus the computational cost of merging cells 
is kept as low as that of non-staggered models. Because the 
momentum equations are only solved on the fully rectangular 
cells that do not enclose the entire huid domain, momentum is 
not conserved in our model as it is in many other models using 
a vertically staggered grid (e.g., Saito et al. 2001; Satoh 2002; 
Klemp et al. 2007). Our model also does not conserve kinetic 
energy, although mass-conservation is guaranteed. 

The cut-cell calculation described in this subsection is only 
performed in the hnest cubes. The ghost cells that are added 
beyond the boundary of each cube allow merging of cells 
between adjacent hnest cubes. To prevent cell merging be¬ 
tween cells that belong to different sizes of cubes, the size 
of the cubes near the terrain surface is readjusted in the grid- 
generation process, if necessary, so that both merged and non- 
merged cut-cells are inside of the hnest region. In the coarser 
cubes, the standard hnite-difference code is used so that the 


computational overhead induced by the use of cut-cell code 
only occurs in the area near the terrain surface. While on cut- 
cells the scheme is locally hrst-order, a global second-order 
accuracy is achieved, as demonstrated numerically in the ap¬ 
pendix. 


E. Parallelization 

The parallelization of the model is straightforward. Load 
balancing is achieved by distributing an equal number of cubes 
at each rehnement level for each processor. To keep the load 
balance in the presence of cut-cells and underground cells at 
the hnest level, we compute solutions for all cells in the hnest 
cubes using the same cut-cell code, and then set the values at 
cells that are located below the terrain surface to zero. The 
overhead required to calculate values of underground cells 
is sufficiently low because cells in cubes that are completely 
bellow the terrain surface are excluded from the computation 
(Figure Ic). For efficient parallelization, a sufficient number of 
cubes would be required to avoid inequality in the distributed 
numbers of cubes among processors. 

The model code is parallelized using OpenMP. The scala¬ 
bility of the model is demonstrated in section IIIB through 
mountain-wave simulations both on a uniform grid and a lo¬ 
cally rehned grid. The data structure of the block-structured 
mesh used in this study can also be suitable for parallel 
computing by Message Passing Interface (MPI) in a dis¬ 
tributed memory system, as demonstrated using BCM by 
Takahashi et al. (2008) and Sakai et al. (2013). 


III. RESULTS 

This section presents the model results from test simulations 
of How passing over a 3D hill. Every test involves an isolated 
hill located at the center of the lower boundary. The model is 
integrated for 1 hour for each test. The following atmospheric 
and boundary conditions are imposed on all the simulations. 
A constant horizontal velocity and Brunt-Vaisala frequency, 
(7 = 10 m s“^ and N = 0.01 s“^, respectively, are initially 
imposed on the entire domain. A sea level potential temper¬ 
ature is specified to be 0 = 300 K. The lower and lateral 
boundary conditions are free-slip and cyclic, respectively. To 
prevent cyclic lateral boundaries from contaminating the sim¬ 
ulated results, a large domain length of 64 km is used in the 
X direction. In the y direction, the domain length is set to 32 
km. The height of the domain is 16 km, and a sponge layer 
(Klemp and Lilly 1978) is placed higher than 10 km to avoid 
rehecting the gravity wave at the rigid top boundary. 


A. Cut-cell model on a uniform mesh 

First, we test the model on a uniform grid of equal spacing 
in order to validate the cut-cell representation of topography 
without using any grid rehnement. The performance of the 
model on a locally rehned grid is discussed in section IIIB. 
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1. Flow over a 3D bell-shaped hill 


This test examines the accuracy of the cut-cell representa¬ 
tion of topography by simulating a flow over a 3D bell-shaped 
hill and comparing the result to that predicted by theory. The 
surface height of the hill is described as, 


h{x,y) 


(1 + jo? -f ’ 


(59) 


where the maximum height = 400 m and the half-width 
0=1 km. Following Locket al. (2012), the general shape of 
the hill is chosen to agree with that analyzed in Smith (1980). 
Here the same grid spacing of 200 m is used in x, y and z 
directions over the entire domain. 

The maximum slope angle of the hill resolved by the model 
grid is 18.1°. The horizontal gradients of the terrain surface, 
calculated by Eqs (17) and (18), are less than 1 at all grid 
columns, so that only the vertical merging of cells are used 
in this test. 

Figure 10 shows the simulated vertical velocity fields over 
the 3D bell-shaped hill. The results successfully demonstrate 
stable and smooth solutions over many integrations. A vertical 
x-z cross section through the center of the hill (Figure 10(a)) 
shows the generation of vertically propagating waves from the 
top of the hill. The phase lines tilt upstream with increasing 
height above the terrain surface, just as they do in the 2D flow 
over a bell-shaped hill with the same height and half-width, 
such as examined in Gallus and Klemp (2000). Unlike the 2D 
case, the amplitude of the flow decays rapidly with height be¬ 
cause of the rapid reduction of the energy density due to the 
three-dimensionality. Figures 10(b) and (c) show the horizon¬ 
tal x-y slices of the vertical velocity field. Above the surface at 
z = 800 m (Figure 10(b)), the fluid rises on the upstream side 
of the hill and descends over the leeward, showing a damped 
oscillation towards downstream. The wavelength of the oscil¬ 
lation is about 7 km, which is comparable to that calculated by 
the linear theory in Smith (1980), 2ttU/N 6.3 km, and also 
agrees well with the numerical result by Lock et al. (2012). At 
the higher altitude of z = 2000 m (Figure 10(c)), the flow pat¬ 
terns on the leeward of the mountain splits and widens to form 
U-shaped regions in good agreement with the theory. 

These results lead to the conclusion that the model repro¬ 
duces the 3D structure of flow over the bell-shaped hill, includ¬ 
ing detailed features as predicted by the theory. In addition, it 
is shown that vertical combinations of cells in our model work 
properly over varying gradients. 


2. Flow over a 3D hemisphere-shaped hill 

The second test demonstrates the ability of the cut-cell 
model to simulate flow over very steep topography, using a 
hemisphere hill of radius r = I km. As in the case with the 
bell-shaped hill, a uniform grid is used in this test with equal 
spacing of 200 m in all three directions. The maximum slope 
angle is resolved at 71° at the foot of the hill. Because the 
horizontal gradients of the terrain surface range from gentle 
around the top of the hill to very steep at the foot, all five di¬ 
rections of cell-merging described in section IIC are involved 
in this test. 



25 30 35 40 45 50 55 
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FIG. 10. The vertical velocity over the 3D bell-shaped hill after 1 
hour of integration: (a) a vertical x-z slice through the center of the 
hill, (b) a horizontal x-y slice at height z = 800 m and (c) a horizontal 
x-y slice at height z = 2000 m. Contour intervals are 0.25 m s“^ in 
(a) and 0.1m s“^ in (b) and (c). Solid and dashed lines indicate 
positive and negative values, respectively. Shaded regions in (b) and 
(c) indicate the region within the half-width from the top of the hill. 


Figure 11(a) shows a vertical x-z cross section of the simu¬ 
lated vertical velocity held. Due to the steepness of the topog¬ 
raphy, much larger amplitudes are calculated in the flow over 
the hemisphere hill compared to the result of Figure 10(a) for 
the bell-shaped hill. Figures 11(b) and (c) show the horizontal 
x-y slices of the vertical velocity held over the hemisphere hill 
for heights of 1200 m and 2400 m, respectively. In the leeward 
side of the hill, damped oscillations towards downstream and 
U-shaped flow patterns widening with height are observed as 
in the case of the bell-shaped hill. Note that contour intervals 
in Figures 11(b) and (c) are two and a half times larger than 
that in Figures 10(b) and (c). 

When the flow past an isolated sphere is symmetric about 
a horizontal plane through the centre, it can be equated to the 
flow past a hemisphere with the free-slip boundary condition 
at the surface (Baines 1995). Numerical results by Hanazaki 
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FIG. 11. Results of flow over the 3D hemisphere hill after 1 hour of 
integration, (a) Vertical x-z slice of the vertical velocity field through 
the center of the hill, (b)-(c) Horizontal x-y slices of the vertical 
velocity field at height 2 = 1200 m and 2 = 2400 m, respectively, (d) 
Vertical x-z slice of the potential temperature field through the center 
of the hill. Contour intervals are 0.25 m s“^ in (a), (b) and (c), and 1 
K in (d). Solid and dashed lines indicate positive and negative values, 
respectively. Shaded regions in (b) and (c) indicate the position of the 
hill. 


(1988) indicate that the flow past a sphere for the case with 
Nr/U = 1 is indeed symmetric about the horizontal plane. 
Referring to the isopycnic lines over the sphere calculated in 
Hanazaki (1988) (see their Figures 3(c) and 11(f)), our model 
faithfully reproduces a smooth flow field close to the terrain 
surface, as illustrated in the result of the potential tempera¬ 
ture on a vertical x-z cross section through the center of the 


hill (Figure 11(d)). Lin et al. (1992) examined experimentally 
the lee wavelength in the flows past a sphere. The observed 
wavelength in the lee waves for the case with Nr/U = 1 was 
slightly less than 2ttU/N « 6.3 km, which is comparable to 
the apparent lee wavelength in our model result of Figure 11 (b) 
(approximately 7 km). 

Though the maximum slope angle of the hemisphere hill is 
almost as large as that of a very steep bell-shaped hill used in 
Locket al. (2012), the reproduced flows show very different 
structures; vertically propagating mountain waves are repro¬ 
duced in this test, whereas a simple flow that climbs over and 
goes down the hill is reproduced in Lock et al. (2012). These 
results can be explained in terms of the parameter Na/U 
indicating the types of disturbances that would be observed 
in each case, according to linear theory (e.g., Satomura et al. 
2003, Smith 1977). When Na/U 3> 1, the theory pre¬ 
dicts that nearly hydrostatically balanced mountain waves will 
appear, while nonhydrostatic mountain waves should be ob¬ 
served when Na/U w 1. Evaluating Na/U by Nr/U = 1 
for the case of the hemisphere hill, the model results in Fig¬ 
ure 11 show the generation of mountain waves as expected 
from the theory. Given the steep nature of the hemisphere hill, 
we may not expect the linear theory to give a good descrip¬ 
tion of the flow. Flowever, referring to the flow past a sphere 
when Nr/U = 1, which is classified by Lin et al. (1992) as 
their flow regime “WV”, our result appears to be compara¬ 
ble to the linear theory as discussed in Baines (1995). When 
Na/U 1, buoyancy forces become less important and the 
motion approaches potential flow, such as observed in the case 
of Locket al. (2012). 

We therefore conclude that our model successfully repro¬ 
duces a smooth and physically reasonable flow over a hill 
which has very steep slopes. At the same time, it is proved 
that the model works properly even when vertically and hori¬ 
zontally merged cells coexist. 


B. Cut-cell model on a locally refined mesh 

Finally we demonstrate the performance of our cut-cell 
model using the block-structured mesh-refinement approach 
described in section IIB. In this test, four simulations of 
flow over the same 3D hemisphere-shaped hill used in section 
III A 2 are performed using grids with different refinement lev¬ 
els. 

The grids used in this test are shown in Figure 12. Each fig¬ 
ure shows the boundaries of cubes that comprise the grid. Top 
panels describe the cube boundaries on a vertical x-z cross 
section through the center of the hill, and bottom panels de¬ 
scribe those on a horizontal x-y cross section at the sea level. 
Here every grid contains 20 x 20 cells in each cube. At the re¬ 
finement level 0 (Eigure 12(a)), the domain is filled with large 
cubes of 400 m spacing in all three directions. At the refine¬ 
ment level 1 (Eigure 12(b)) and 2 (Eigure 12(c)), cubes are 
divided into finer cubes around the hill, where the medium- 
and small-sized cubes have spacings of 200 m and 100 m in 
all three directions, respectively. The results using the grids 
with each refinement level are compared to the result using 
a uniformly fine mesh with equal spacing of 100 m (Eigure 
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FIG. 12. Computational grid used in the test of flow over the 3D hemisphere hill with refinement levels (a) 0, (b)l, (c) 2 and with (d) a 
uniformly fine mesh. The top and bottom panels in each figure describe the vertical and horizontal structures of the grid. Solid lines are the 
cube boundaries and each cube has 20 x 20 cells. The grid spacing in large, medium and small cubes is 400 m, 200 m, and 100 m, respectively. 
Shaded region describes the topography. 


12(d)), which is considered as a reference. The total numbers 
of cubes are 64, 148, 232 and 4096 in Figures 12(a), (b), (c) 
and (d), respectively. 

Figures 13(a) to (d) respectively show vertical x-z cross sec¬ 
tions of the simulated vertical velocity fields using the grids 
with refinement levels 0, 1, 2 and using the uniformly fine 
mesh. Due to the lack of resolution, the amplitude of the ver¬ 
tical velocity using the grid with the refinement level 0 (Fig¬ 
ure 13(a)) decays faster with height compared to that using 
the uniformly fine mesh (Figure 13(d)). In the leeward side 
of the hill, the propagating waves towards further downstream 
captured in Figure 13(d) are missing in Figure 13(a). The am¬ 
plitude is successfully increased when the refined meshes are 
implemented around the hill (Figures 13(b) and (c)). Though 
in the region of the coarse resolution the simulated amplitude 


is smaller than that in the reference solution of Figure 13(d), 
both results using the grid with the refinement level 1 (Figure 
13(b)) and the refinement level 2 (Figure 13(c)) successfully 
reproduce the wave propagation towards downstream. In ad¬ 
dition, no visible distortions of mountain wave patterns are 
found in these results across the fine-coarse cube boundaries. 

To compare the results of Figures 13(b) and (c) in detail, 
the differences with respect to the reference solution of Figure 
13(d) are shown in Figures 13(e) and (f), respectively. Note 
that the contour interval in Figures 13(e) and (f) is two fifths 
of that in Figures 13(a)-(d). In Figure 13(e), some errors 
are found in both the regions at the refinement levels 0 and 
1 where the resolution is coarser than that of the reference so¬ 
lution. These errors are successfully reduced when implement¬ 
ing the refined mesh at level 2 around the hill in Figure 13(f). 
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FIG. 13. Vertical x-z slices of the vertical velocity field through the center of the 3D hemisphere hill in the mesh refinement test case. Results 
are reproduced on the grid with refinement levels (a) 0, (b)l, (c) 2 and with (d) a uniformly fine mesh. Figure (e) shows the difference between 
the vertical velocity fields in (b) and (d), and (f) shows the difference between (c) and (d). Contour intervals are 0.25 m s“^ in (a)-(d), and 0.1 
m s“^ in (e) and (f). 


In particular, it is shown that refining the mesh around the hill 
not only improves the result in the refined region at level 2 but 
also improves the result in the coarser regions at levels 1 and 
0. These results indicate that the benefit due to the higher res¬ 
olution within the refined cubes exceeds the errors due to the 
resolution changes at fine-coarse cube boundaries. 

Figure 14 displays horizontal x-y slices of the simulated 
vertical velocity fields. Here Figures 14(a) and (b) show the 
vertical velocity fields at z = 1200 m using the grid with the 
refinement level 2 (Figure 12(c)) and using the uniformly fine 
mesh (Figure 12(d)), respectively. Figure 14(c) shows the dif¬ 
ference between the velocity values in Figures 14(a) and (b), 
where the contour interval in Figure 14(c) is two fifths of that 
used in Figures 14(a) and (b). As in the vertical slices of the 
results, the model reproduces very accurate flow in the region 


at the refinement level 2 compared to the reference solution 
using the uniformly fine mesh. As the mesh resolution gets 
coarser towards the leeward side of the hill, the amplitude de¬ 
cays faster in Figure 14(a) compared to that in the reference 
solution of Figure 14(b), therefore the errors become larger 
as shown in Figure 14(c). Because there are no visible distor¬ 
tion nor large errors found at the fine-coarse mesh boundaries, 
these errors are considered to be mainly due to the reduction of 
the mesh resolution. A similar behavior of the model is found 
at z = 2400 m (Figures 14 (d) to (f)). Here Figures 14(d) and 
(e) show the vertical velocity fields at z = 2400 m using the 
grids of Figures 12(c) and (d), respectively, and Figure 14(f) 
shows the difference between them. These results confirm the 
model’s ability to effectively increase the resolution around 
the hill without introducing large errors associated with the 
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(a) refinement level 2 


(d) refinement level 2 



(b) uniformly fine mesh (e) uniformly fine mesh 





X (km) X (km) 


FIG. 14. Horizontal x-y slices of the vertical velocity field over the hemisphere hill at height « = 1200 m reproduced on the grid with (a) the 
refinement level 2 and (b) a uniformly fine mesh. Figure (c) shows the difference between the vertical velocity in (a) and that in (b). Solid and 
dashed lines indicate positive and negative values, respectively. Contour intervals are 0.25 m s“^ in (a) and (b), and 0.1 m s“^ in (c). Gray lines 
in (a) and (c) describe the cube boundaries with the refinement level 2. Shaded region describe the topography, (d)-(f) same as with (a)-(c) 
except 2 = 2400 m. 


change of mesh resolution. 

Finally the speed-up ratio is estimated in order to verify the 
scalability of the model both on a locally refined grid and on 
a uniform grid (Figure 15). The parallelization of the model 
code is described in section HE. Here the computations are 
performed on a Linux PC which has 16-core AMD Opteron 
CPUs. The parallel efficiency is about 80% and 70% at 8 and 
16 processors, respectively, in both tests using the grid with 
the refinement level 2 (Figure 12(c)) and using the uniformly 
fine mesh (Figure 12(d)). This result indicates that the mesh 
refinement and the subcycling time integration do not affect 
the parallel efficiency of the overall computation in our model. 

These results lead to the conclusion that the model repro¬ 
duces a flow result on a locally refined grid with both sufficient 


accuracy and computational efficiency, thereby demonstrating 
the advantage of the model’s mesh-refinement technique in 
combination with the cut-cells for high-resolution atmospheric 
simulations over 3D topography. 


IV. CONCLUSIONS 

To achieve high-resolution and highly precise simulations 
over 3D topography, a new 3D nonhydrostatic atmospheric 
model was developed using the Cartesian cut-cell method. To 
avoid severe restriction on time steps while maintaining the 
sharp representation of terrain surface, merged cells are used 
in both the horizontal and vertical directions. In addition, a 


















































































17 



FIG. 15. Comparison of speed-up ratio with increasing the number 
of threads. 


block-structured mesh-refinement approach was implemented 
to allow high resolution at the boundary layers with computa¬ 
tional efficiency. 

The accuracy of the model was confirmed by simulating 
flow over a 3D bell-shaped hill and comparing the result with 
the analytical solution given by the linear theory. The abil¬ 
ity of the model to simulate flow over topography including a 
very steep slope was demonstrated using a hemisphere-shaped 
hill. The model successfully reproduced smooth and physi¬ 
cally reasonable flow including detailed features which agree 
with numerical and experimental results. Finally, a test of the 
mesh-refinement approach in combination with the cut-cell 
representation of topography was performed using the same 
hemisphere-shaped hill. The model reproduced smooth moun¬ 
tain waves propagating over varying grid resolution without 
introducing large errors at fine-coarse cube boundaries. In ad¬ 
dition, the model showed a good parallel efficiency on a lo¬ 
cally refined grid around the hill. 

A remaining issue with the use of Cartesian cut-cell grids is 
the incorporation of physical process models and other compo¬ 
nents, such as a land surface model. As many of the physical 
parameterization schemes are based on a single column model, 
the couplings with the block-structured grid and with horizon¬ 
tally merged cells may not be straightforward. Further work 
is intended to solve these issues with the aim of simulating 
real-world atmospheric problems. 


ACKNOWLEDGMENTS 



1.0E+02 . , , 1.0E+03 

A(m) 

FIG. 16. Variation of global error in u and w with grid intervals for 
flow over the cosine-shaped hill. Thin solid and dashed lines indi¬ 
cate Li and L 2 norms of the errors, respectively, and thick solid line 
corresponds to second-order accurate convergence. 


This work was supported by the JSPS Postdoctoral Fellow¬ 
ship for Research Abroad. Part of this work was supported 
by a Grant-in-Aid for JSPS Bilateral Joint Research Project. 
Some figures were drawn by the GFD-DENNOU Library. 


Appendix A: Estimation of Globai Accuracy 


In this section, we numerically demonstrate that our cut-cell 
method is globally second-order accurate for a flow over 3D 
terrain. Here we performed four simulations of flow over a 3D 
cosine-shaped hill using a uniform grid with different grid in¬ 
tervals: A = 400 m, 250 m, 200 m, and 100 m. The same grid 
interval is used in x, y and z directions for each test, and the 
same atmospheric and boundary conditions as used in section 
III are imposed. The shape of the hill is described by 


( 


1 0 



if r < a 
if r A a, 


(Al) 
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where 

r = ^{x- Xc)2 + {y - VcY (A2) 

denotes the distance between a position {x, y) and the center 
of the lower boundary (xc, yc)'- in this test, Xc = 32 km and yc 
= 16 km. The maximum height /im and the half-width a are 
set to 400 m and 2 km, respectively. 

Eollowing Yamazaki and Satomura (2010), both the Li and 
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L 2 norms of the global errors are computed as, 


= 


^L2 


N^NyN, 


numerical „i,exact 


1 




NM.N., 


x^yy^^z j_ 


y (V': 


-i’T 

_ ^exacty^ 


numerical ^i,exact\2 


7=1 


(A3) 

1/2 

5 

(A4) 


respectively, by assuming the solution with the finest solution 
of Z\ = 100 m as the exact solution. Here N^, Ny and 
indicate the grid numbers above the top of the hill in x, y and 
z directions, respectively. Figure 16 shows a log-log plot of 
the errors in velocity components u and w versus grid intervals. 
Also shown is a line with a slope of 2 which corresponds to the 
second-order accurate convergence. The plot clearly shows 
that the global error in our computed solution decreases in a 
manner consistent with a second-order accurate scheme. This 
test therefore shows that our method produces results which 
are consistent with a method of global second-order accuracy. 
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